function [NLI_sim_data,NLI_sim_data_3d] = NLI_sim_funct(vp,fp,col_long,data)

% Unpack parameters of two-step logit specification:

% (1) logit for prob(I>0) = 1 - (exp(.)/(1+exp(.)))
mu1 = vp(1); % intercept
mu2 = vp(2); % mom age, positive
mu3 = vp(3); % mom agesq
mu4 = vp(4); % mom educ, positive
mu5 = vp(5); % dad age, positive
mu6 = vp(6); % dad agesq
mu7 = vp(7); % dad educ, positive

% (2) I = exp(intercept + ...)
mu8 = vp(8); % intercept
mu9 = vp(9); % mom age, positive
mu10 = vp(10); % mom agesq
mu11 = vp(11); % mom educ
mu12 = vp(12); % dad age, positive
mu13 = vp(13); % dad agesq
mu14 = vp(14); % dad educ
sigma_NLI = vp(15); % sigma

NLI_sim_data = nan(fp.N*(fp.M+1)*fp.R,fp.tot_cols);
NLI_sim_data_3d = nan(fp.M+1,fp.N,fp.R); % useful for calling whenever we use backward induction

sim_data_index = 0;
for r = 1:fp.R
    for q = 1:fp.N %this should be same as fp.N
        sim_id = r*1000 + q; % don't use i instead of q! We need a unique identifier for every simulated household
        % (Note: if we would have more than 1000 hhs, use r*10000 + q)

        % mother/child i's data
        data_here = data(1+(fp.M+1)*(q-1):(fp.M+1)*q, : ); % note: this could be bootstrapped data

        % Load exogenous household characteristics
        if fp.compstat_SES_homwages == 1
            % In comparative statics exercise, shut down heterogeneity in wages and NLI by parental education level
            mom_educ = fp.s1_median;
            dad_educ = fp.s2_median;
        elseif fp.compstat_SES_homwages == 0
            mom_educ = data_here(1,col_long.mom_educ); % constant over time
            dad_educ = data_here(1,col_long.dad_educ);
        end
        if fp.compstat_age_homwages == 1
            % In comparative statics exercise, shut down heterogeneity in wages and NLI by parental age
            mom_age = fp.age1_median * ones(fp.M+1,1);
            dad_age = fp.age2_median * ones(fp.M+1,1);
        elseif fp.compstat_age_homwages == 0
            mom_age = data_here(:,col_long.mom_age); % vector of size (fp.M + 1) * 1
            dad_age = data_here(:,col_long.dad_age); % vector of size (fp.M + 1) * 1
        end

        for t = 0:fp.M
            % Note: t corresponds to the child's age, so we pass this along in the stackelberg solver
            sim_data_index = sim_data_index + 1; % we use this index to fill out the simulated data matrix

            % We don't simulate for years prior to 1997, but fill out the id  and exogenous
            % variables (of the simulated household) nonetheless (for the moment function)
            NLI_sim_data(sim_data_index,col_long.id) = sim_id;
            NLI_sim_data(sim_data_index,col_long.mom_age) = mom_age(t+1); % this is a (fp.M+1)*1 vector
            NLI_sim_data(sim_data_index,col_long.dad_age) = dad_age(t+1); % this is a (fp.M+1)*1 vector
            NLI_sim_data(sim_data_index,col_long.mom_educ) = mom_educ;
            NLI_sim_data(sim_data_index,col_long.dad_educ) = dad_educ;
            unif_inc = fp.NLI_unif_draws(q,r,t+1,1); % uniform "logit" draw
            epsilon_inc = fp.wage_draws(q,r,t+1,3); % normal epsilon draw

            % baseline model: two-step with logit
            tmp1 = exp(mu1 + mu2*mom_age(t+1) + mu3*(mom_age(t+1)^2) + mu4*mom_educ + mu5*dad_age(t+1) + mu6*(dad_age(t+1)^2) + mu7*dad_educ);
            tmp2 = (tmp1/(1+tmp1)); % probability of success, aka I>0
            % now compare this to our random number draw:
            tmp3 = (unif_inc <= tmp2); % hence, tmp3 = 1 with probability tmp2, which is what we want (and unif_inc is kept fixed by a seed 'state')
            if tmp3 == 0
                % in this case, I = 0;
                tmp = 0;
            elseif tmp3 == 1
                % in this case, I has to be strictly positive
                tmp = exp(mu8 + mu9*mom_age(t+1) + mu10*(mom_age(t+1)^2) + mu11*mom_educ + mu12*dad_age(t+1) + mu13*(dad_age(t+1)^2) + mu14*dad_educ + sigma_NLI * epsilon_inc);
            end

            if tmp>1000
                tmp = 1000; % censor these observations (note that we cannot have a NaN observation in the estimation!)
            end
            % save our generated NLI value (in 2 different formats for later use)
            NLI_sim_data(sim_data_index,col_long.nonlabor_income) = tmp;
            NLI_sim_data_3d(t+1,q,r) = tmp;
        end
    end
end
end